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Abstract 

When neurons fire action potentials, dissipation of free energy is usually not directly considered, because the change in free 
energy is often negligible compared to the immense reservoir stored in neural transmembrane ion gradients and the long- 
term energy requirements are met through chemical energy, i.e., metabolism. However, these gradients can temporarily 
nearly vanish in neurological diseases, such as migraine and stroke, and in traumatic brain injury from concussions to severe 
injuries. We study biophysical neuron models based on the Hodgkin-Huxley (HH) formalism extended to include time- 
dependent ion concentrations inside and outside the cell and metabolic energy-driven pumps. We reveal the basic 
mechanism of a state of free energy-starvation (FES) with bifurcation analyses showing that ion dynamics is for a large 
range of pump rates bistable without contact to an ion bath. This is interpreted as a threshold reduction of a new 
fundamental mechanism of ionic excitability that causes a long-lasting but transient FES as observed in pathological states. 
We can in particular conclude that a coupling of extracellular ion concentrations to a large glial-vascular bath can take a role 
as an inhibitory mechanism crucial in ion homeostasis, while the Na+/K+ pumps alone are insufficient to recover from 
FES. Our results provide the missing link between the HH formalism and activator-inhibitor models that have been 
successfully used for modeling migraine phenotypes, and therefore will allow us to validate the hypothesis that migraine 
symptoms are explained by disturbed function in ion channel subunits, Na+ZK"*^ pumps, and other proteins that regulate 
ion homeostasis. 
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Introduction 

The Hodgkin-Huxley (HH) model is one of the most successful 
models in mathematical biology [1]. This formalism, i.e., a HH- 
type model, describes voltage changes across cell membranes that 
result in excitability. Not only neurons are excitable cells, also 
myocytes, pancreatic jS-ceUs, and even a plant cell (Chara 
corallina) exhibit excitable dynamics [2-4]. The dynamic range 
of phenomena includes single action potentials (spikes), periodic 
spiking, and bursting (slow modulation of spiking). For example, in 
pancreatic /?-cells bursting is induced by a calcium current [4,5] . 
A more complete treatment of this phenomenon, however, also 
requires the inclusion of Na+/K^ pumps [6]. The dynamics of 
ion pumps and ion concentrations is also crucial for cardiac 
alternans (periodic beat-to-beat variations) and higher-order 
rhythms in the ischemic ventricular muscle [7-9]. 

In the literature such augmented HH-type models are also 
called second-generation HH models [10]. In the context of 
certain pathologies of the brain, whose fundamental dynamic 
structure we study here, we prefer the simpler name 'ion-based' 
models. This indicates that ion concentrations are major 
dynamical, that is, time-dependent variables. Their dynamical 
role in neuron models goes beyond merely modulating spiking 
activity. Ion dynamics can lead to a completely new type of ionic 
excitability and bistability, that is, the phenomena of so-called 



'spreading depolarizations' and 'anoxic depolarization', respec- 
tively. (Spreading depolarizations are also called 'spreading 
depression' and we will use both names interchangeably in this 
paper.) These depolarized states of neurons are related to 
migraine, stroke, brain injury, and brain death, that is, to 
pathologies of the brain in which a transient or permanent 
break-down of the transmembrane potential occurs [11,12]. 
Another even more characteristic property of this 'twUight state 
close to death' [13] are the nearly completely flat transmembrane 
ion gradients. The almost complete break-down of both 
membrane potential and — due to reduced ion gradients — Nernst 
potentials together cause a nearly complete release of the Gibbs 
free energy, that is, the thermodynamic potential that measures 
the energy available to the neurons for normal functioning. We 
hence refer to this state as a state of free energy-starvation (FES). 
We want to stress that such phenomena require the broader 
thermodynamical perspective, because it goes beyond the HH 
description in terms equivalent electrical circuits in membrane 
physiology (see discussion). 

The object of this study is to clarify quantitatively the detailed 
ion-based mechanisms, in particular the time-dependent poten- 
tials, leading to this condition. In fact, early ion-based models 
have been introduced in a different context to describe excitable 
myocytes and pancreatic /?-cells with variable ion concentrations 
[14-16]. Neuronal ion-based models have been used to study 
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Author Summary 

Theoretical neuroscience connplennents experinnental and 
clinical neuroscience. Simulations and analytical insights 
help to interpret data and guide our principal under- 
standing of the nervous systems in both health and 
disease. The Hodgkin — Huxley-formulation of action 
potentials is certainly one of the most successful models 
in mathematical biology. It describes an essential part of 
cell-to-cell communication in the brain. This model was 
in various ways extended to also describe when the 
brain's normal performance fails, such as in migraine 
hallucinations and acute stroke. However, the fundamen- 
tal mechanism of these extensions remained poorly 
understood. We study the structure of biophysical neuron 
models that starve from their 'free' energy, that is, the 
energy that can directly be converted to do work. 
Although neurons still have access to chemical energy, 
which needs to be converted by the metabolism to 
obtain free energy, their free energy-starvation can be 
more stable than expected, explaining pathological 
conditions in migraine and stroke. 



spreading depolarizations (SD) [17-22] and anoxic depolariza- 
tions [21]. In these phenomenological studies the types of ion 
dynamics related to the pathologies have been reproduced, but 
not investigated in a bifurcation analysis. Hence the fundamental 
phase space structure of these high-dimensional models that 
underlies the ionic excitability characterisitic of SD remains 
poorly understood. Furthermore, neuronal ion-based models 
have been used to study seizure activity [23,24] and spontaneous 
spiking patterns in myelinated axons with injury-like membrane 
damaging conditions (e.g., caused by concussions) [25,26]. In 
tiiese models, the phase space structure was investigated, 
however, only with respect to the modulating effect of ion 
concentrations on the fast spiking dynamics (seizure activity, 
injuries), and with respect to spiking node-to-node transmission 
fidelity (myelinated axons). 

In this paper we present bifurcation analyses of several minimal 
biophysical ion-based models that reveal bistability of extremely 
different ion configurations — physiological conditions vs. free 
energy-starvation — for a large range of pump rates. In related 
models certain bistabUities have been explored before. For 
example, Frohlich et al. [27-29] found coexistence of quiescence 
and bursting for certain fixed extracellular potassium concentra- 
tions and also bistability of a physiological and a strongly 
depolarized membrane state in a slow— fast analysis of calcium 
gated channels. Bistability of similar fixed points has also been 
found for the variation of extracellular potassium [30] or, similarly, 
the potassium Nernst potential [31]. Also the effect of pump 
strength variation has been explored under fixed FES conditions 
[32]. In this paper, however, we do not treat slow variables as 
parameters and show bistability of fast dynamics, but instead we 
address the stability of ion concentrations themselves, which are 
subject to extremely slow dynamics. This allows us to find 
bistability of extremely different ion distributions, a feature that 
distinguishes these two states from the polarized and depolarized 
states studied in the afore mentioned work. A study that also had 
significantiy different ion distributions was done by Cressman et al 
[33], however, the seizure-like phenomena discussed in their work 
are quite different — though clinically related — from those pre- 
sented in this paper. 

Because of the occurrence of ion state bistability we conjecture 
that our model describes a threshold reduction of a mechanism 



that leads to ionic excitability in form of spreading depolarizations. 
In other words, we conclude that an important inhibitory 
mechanism to describe ion homeostasis such as glial buffering or 
diffusive regulation of extracellular ion concentrations plays a 
crucial role in ion homeostasis and the Na^ /K^ pumps alone are 
insufficient to recover from free energy-starved states. We show 
that when the extracellular K"'" concentration is regulated by 
linearly coupling it to an infinite bath, the bistable system changes 
to an excitable system, which we call ionic excitability. The effect 
of turning off ghal buffering and diffusion has been discussed in 
more detailed ion-based models [27,29] before, but has not been 
related to the fundamental phase space structure of the system. 
Our conclusions have been validated by demonstrating the 
robustness of the results in a large variety of minimal ion-based 
models, which all consistendy show this insufficiency of Na"'" /K"'" 
pumps, and also in a very detailed membrane model that has been 
used intensively for computational studies of spreading depolar- 
izations and seizure-like activity [17,34]. 

Model 

Hodgkin-Huxley (HH) nnodel and reductions 

A simple ion-based neuron model can be obtained as a natural 
extension of the Hodgkin-Huxley (HH) model [1]. We list the 
basic equations of HH that we used for the sake of completeness, 
and also comment on two often used model reductions of which 
one must be modified for our study. Furthermore leak currents are 
specified, which is necessary for the extension towards ion-based 
modeling. 

In the HH model, single neuron dynamics is described in terms 
of an electrically active membrane carrying an electric potential V, 
and the three gating variables n, m and h that render the system 
excitable. Ion species included are sodium, potassium, and an 
unspecified ion carrying a leak current, which can be attributed to 
chloride in our extended model. The rate equations read [1] : 



d_F 



1 



-(I 



Na + 



-Ik+ +IcI lapp), 



dn 
dt 



dh haj—h 
dt T/, 

dm nias—in 
dt T„ 



(2) 



(3) 



(4) 



The top equation is simply KirchhofFs current law for a 
membrane with capacitance Cm and membrane potential V. 
lapp is an externally applied current that may, for example, 
initiate voltage spikes. The gating variables n, h, and m are the 
potassium activator, sodium inactivator, and sodium activator, 
respectively. Their dynamics is defined by their voltage- 
dependent asymptotic values x^o and relaxation times Tx {x = n, 
in, h). These are given by 



and 
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T,. = — — for x = n, in and h. 

■ <^(=(, + ;8,) 

Here 0 is a common timescale parameter, and the Hodgkin- 
Huxley exponential functions are 



0.1(K + 30) 



l-exp(-(F + 30)/10)' 



ft„=4exp(-(K + 55)/lf 



0.01(K+34) 
l-exp(-(K + 34)/10)' 



i?„ = 0. 1 25 exp( - ( F + 44)/80), 



ai, = 0.07 exp( - ( F + 44) /20)), 



1 

l+exp(-0.1(F+14))' 



(5) 



(6) 



(7) 



(8) 



(9) 



(10) 



The individual ion currents read 



lNa+ =(g'Na+^Na"^^h)iV-ENa), 



la- =g'aiV-Ea), 



(11) 



(12) 



(13) 



with g;* denoting leak and gated conductances. In fact, Hodgkin 
and Huxley set up their model with an unspecified leak current 
and non-leaking sodium and potassium channels. As long as ion 
dynamics is not considered this is mathematically equivalent to 
specifying the leak current as being partially sodium, potassium 
and chloride, but it is physically inconsistent because the reversal 
potentials for the ions differ. In an ion-based approach, however, 
the main task of the ion pumps under physiological conditions is to 
compensate for sodium and potassium leak currents (see next 
section) while gated currents are extremely small in the 
equilibrium. So at this point leak currents for all ion species are 
important. 

The Nernst potentials Eig„ are given in terms of the ion 
concentrations [ion] in the intracellular space (ICS) and the 
extracellular space (ECS) denoted by subscripts / and e, 
respectively: 



26.64 



H[ion\J[ion\.), 



(14) 



for ion = K, Na, and CI and ^,o« is the ion valence. All model 
parameters are listed in Table 1. The units chosen are those 
typically used and appropriate for the order of magnitude of the 



Table 1. Parameters for Hodgkln-Huxley model. 





Name 


Value & unit 


Description 




1 i-LF/cm- 


membrane capacitance 




3/msec 


gating timescale parameter 




0.0175 mS/cm^ 


sodium leak conductance 


s%„ 


100 mS/cm^ 


max. gated sodium conductance 


s'k 


0.05 mS/cm^ 


potassium leak conductance 


sic 


40 mS/cm^ 


max. gated potassium conductance 


S'a 


0.05 mS/cm^ 


chloride leak conductance 


Noi 


27 niMol// 


ECS sodium concentration 


Na, 


120 mMol// 


ICS sodium concentration 


Ki 


130.99 mMol// 


ECS potassium concentration 


Ke 


4 mMol// 


ICS potassium concentration 


Cli 


9.66 mMol// 


ECS chloride concentration 


Ch 


124 mMol// 


ICS chloride concentration 




39.74 mV 


sodium Nernst potential 


Ek 


-92.94 mV 


potassium Nernst potential 


Ea 


-68 mV 


chloride Nernst potential 



doi:1 0.1 371 /journal.pcbi.l 003551 .tool 

respective quantities. Time is measured in msec, potentials in mV, 
and ion concentrations in mMol//. The units for conductance 
densities imply that ionic and pump current densities are in 
(tA/cm^. For better readability we omit the square brackets on the 
ion concentrations and simply write Kij^, Naijg, and C/,yp. 

For lapp = 0 this model is monostable with an equilibrium at 
K = — 68 mV. Note that EffaJ^V and Ekj^V imply that under 
equilibrium conditions neither /jVff+ nor Ix+ vanish, but only their 
sum does. Sufficiently strong current pulses can — depending on 
their duration — initiate single voltage spikes or spike trains. 
Constant apphed currents can drive the system to a regime of 
stationary oscillations. The minimal current required for this is 
usually called rheobase current. 

The HH model can be reduced to two dynamical variables in a 
way that preserves these dynamical features. One common 
simplification [35] is to eliminate the fastest gating variable m 
adiabaticaUy and set 



--m^(V). 



(15) 



Second, there is an approximate functional relation between h and 
n that is usually realized as a linear fit [36] . The ion— based model 
presented in this article, however, contains a stable fixed point 
with large n, and a linear best fit would then lead to a negative h. 
Therefore we will use the following sigmoidal fit to make sure h is 
non-negative: 



h = h„g(n)=\- 



I 



l+exp(-6.5(«-0.35)) 



(16) 



After this reduction the remaining dynamical variables are V and 
n. 

Minimal ion-based model 

While in the original HH model ion concentrations are model 
parameters, in ion-based modeling intra- and extracellular ion 
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concentrations become dynamical variables, which causes the 
Nernst potentials to be dynamic. The model defined by the rate 
eqs. (1), (2) and contraint eqs. (15), (16) can straightforwardly be 
extended to make ion concentrations dynamic since currents 
induce ion fluxes. However, under those equilibrium conditions 
found in HH neither =0 nor If^n^ =0. Hence we need to 
include ion pumps [15] to make sure that the rate of change in ion 
concentration inside the cell (/) and extracellular [e) can vanish in 
the resting state (Najj^ = Kjf^, = Cljj^ = 0). 

The rate equations for ion concentrations in the intracellular 
space (ICS) are then 



dNoi 



CO, 



da 

dt 



CO, 



-la- . 



(17) 



(18) 



(19) 



The factor y converts currents to ion fluxes and depends on the 
membrane surface Ai„ and Faraday's constant F: 



Am 



(20) 



Dividing the ion fluxes by the ICS volume co, gives the change 
rates for the ICS ion concentrations. The pump current Ip 
represents the ATP-driven exchange of ICS sodium with 
potassium from the extracellular space (ECS) at a 3/2-ratio. It 
increases with the ICS sodium and the ECS potassium concen- 
tration. Chloride is not pumped. We are using the pump model 
from [23,24]: 



Ip{Na„K,) = p{\ 



-exp( )) 



(l+exp(5.5-i:,))- 



(21) 



where p is the maximum pump current. As a consequence of mass 
conservation ion concentrations in the ECS can be computed from 
those in the ICS [33]: 



Na, = Naf^+'^(Naf^-Na,), 
co„ 



We 



CL 



COp 



-a). 



(24) 



(25) 



(26) 



Dynamic ion concentration imply that the Nernst potentials in eqs. 
(11)-(13) are now dynamic (see eq. (14)). The additional 
parameters of the ion-based model are listed in Table 2. The 
morphological parameters and co,- are taken from [17]. In 
cortical ion-based models, the extracellular volume fraction 
f = 0Je/(0Je + 0Ji) ranges from 13% in [17] to 33% in [21]. In 
experimental studies, / is about 20%, a value that can increase, for 
example, in focal cortical dysplasias type II, a frequent cause of 
intractable epilepsy, to 27% [37] or during sleep to 32% (the latter 
only, if we transfer the increase observed in mouse data to human) 
[38]. It is important to note that in experimental studies, the 
extracellular volume fraction refers to the fraction with respect to 
the whole tissue, which includes also the glial syncytium. Assuming 
equally sized neuronal and glial volume fractions of 40% each, an 
experimentally measured value of 20% would in our model, which 
does not directly include the volume of the glial syncytium, 
correspond to / = 0.33 or 33%. We choose an intermediate value 
of 25% for /, but address the influence of the volume ratio in Sec. 
Results. We prefer to give these morphological parameters in the 
commonly used units which are appropriate to their order of 
magnitude rather than unifying all parameters, e.g. the cell volume 
is given in |tm' instead of / which ion concentrations are related 
to. Consequently y from Table 2 must be multiplied by a factor of 
10 to correctly convert currents to change rates for ion 
concentration in the given units. Because of the extremely small 
value of)' the membrane dynamics, i.e., the dynamics of V and n, 
is five orders of magnitude faster than the ion dynamics. 

A consequence of this large timescale separation is that the 
system will attain a Donnan equilibrium when the pumps break 
down. The Donnan equilibrium is a thermodynamic equilibrium 
state (not to be confused with merely a frxed point, though it is one) 
that is reached for ion exchange across a semipermeable 
membrane. Since we have not explicitly included large imperme- 
able anions inside the cell, this is at first surprising. For no applied 
currents and Ip = 0, the ion rate equations imply that an 
equilibrium requires all ion currents to vanish. Since conductances 



JO) 



, f (0) ■ \ 
H (ion) —lorii}, 

CO,, 



(22) 



with the ECS volume a>g. Superscript zero indicates initial values. 
Since all types of transmembrane currents, i.e., also the pumps, 
must be included in eq. (1) for the membrane potential, we have to 
add the net pump current Ip: 



dV 



1 



(I 



Na + 



-la- +Ip — Iupp)- 



(23) 



The rate equations for the ion-based model are thus given by eqs. 
(2), (17)-(19), (23). These rate equations are complemented by the 
gating constraints eqs. (15), (16) and the mass conservation 
constraints 



Table 2. Model parameters for Ion-based model only. 





Name 


Value & unit 


Description 


CO,- 


2,160 nm' 


volume of ICS 




720 nm' 


volume of ECS 


F 


96485 C/Mol 


Faraday's constant 


A,„ 


922 nm- 


membrane surface 


y 


Lim-Mol 
9.556e 3 ^ 


conversion factor 


p 


5.25 |.iA/cm- 


max. pump current 



doi:l 0.1 371 /journal.pcbi.l 003551 .t002 
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Figure 1. Upper panels: Membrane and Nernst potentials, lower panel: Ion concentrations vs time, (a) Response of the model to a 
0.5 sec long sodium current pulse with amplitude 150 nA/cm^ (marked by the black star). The pulse causes voltage spiking that stops in a strongly 
depolarized state (see blow-up inset). The membrane potential V takes a final value of about — 25mV (upper panel). The ion gradients, i.e., the 
differences between intra- and extracellular ion concentrations, reduce drastically during the stimulation and slowly adjust to a new fixed point after 
a couple of hundreds of seconds (lower panel), (b) Switching off the ion pump for 20 sec (indicated by the light grey interval) causes similar dynamics. 
The membrane depolarization and dissipation of ion gradients is a bit slower than for (a). After the pump is switched on again the system attains the 
same fixed point as in (a). 
doi:10.1371/journal.pcbi.1003551.g001 



are strictly positive it follows that all Nernst potentials and the 
membrane potential must be equal. Ion concentrations will then 
adjust accordingly. However, eqs. (17)-(19) and (23) imply the 
following constraint on the ICS charge concentration Q,: 

AQi : = ANoi + AKi - AC/,- = — A V, (27) 

OJ, 

where A denotes the difference between the initial and final value 
of a variable. Since y is very small, changes in ion concentrations 
must practically satisfy electroneutrality. This condition together 
with the equality of all Nernst potentials defines the Donnan 
equilibrium, so we see that it is contained in our model as the limit 
case with no pumps and no apphed currents. It should be noted 
that this observation provides a necessary condition for the 
correctness of biophysical models. 

In this extension of the HH model the ion dynamics makes 
Nernst potentials time-dependent. The simultaneous effect of a 
diffusive and an electrical force acting on a solution of ions is 
described more accurately by the Goldman-Hodgkin-Katz 
(GHK) equation though. Nevertheless we prefer Nernst 
currents, because this formulation allows us to use well- 
established conductance parameters so that the model is 



completely defined by empirically estimated parameters. In 
Sec. Results we will see how GHK currents can be modelled 
and that the qualitative dynamical behaviour of the system is 
not affected. 

Results 

Phase space analysis of ion-based model 

In the ion — based model introduced above current pulses can 
still initiate voltage spikes (not shown). However, extremely strong 
pulses, in fact comparable to those used in [17] to trigger 
spreading depolarizations, can drive the system away from the 
physiological equilibrium to a second stable frxed point that is 
strongly depolarized (see Fig. 1(a)). This is a new dynamical 
feature. The depolarized state can also be reached when the ion 
pumps are temporarily switched off (see Fig. 1(b)). Apart from the 
depolarization this state is characterized by almost vanishing ion 
gradients. This free energy-starvation (FES) is reminiscent of the 
Donnan equilibrium. Extracellular potassium is increased from 4 
to more than 40mMol/I while the extraceUxilar sodium concen- 
tration is reduced from 120 to less than 30mMol/l. The gated ion 
channels are mostly open (potassium activation n is 60%), and it is 
no longer possible to initiate voltage spikes. In this section we wiU 
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present a phase space analysis of the model and derive conditions 
for the observed bistability between a physiological equilibrium 
and a state of FES. 

Note that the transition from the physiological state to FES 
happens via ion accumulation due to spiking, and we wiU see in 
Sec. Results that indeed the membrane ability to spike is a 
necessary condition for the bistability. Similar processes of ion 
accumulation were regarded as unphysiological in modelling of 
cardiac cells [8], but are familiar in cortical neurons where ion 
accumulation is central to seizure-like activity [24,33] and 
spreading depression [17] (SD). In fact, we will briefly demonstrate 
how the bistability relates to local SD dynamics. 

Symmetry of the ion-based model. Prior to a bifurcation 
analysis we need to discuss a conservation law (symmetry) of eqs. 
(17)-(19), (23). The direct extension of a membrane model to 
include ion dynamics as presented above naturally leads to a linear 
dependence of dynamical variables. In our case this is reflected by 
the following relation 



C, 



dV_ 

' d7 



V 



dNoi dKi 



dt 



+ 



dt 



dCh 

~dr 



(28) 



for lapp = 0. As a consequence the determinant of the Jacobian is 
always zero and the system is nowhere hyperbolic. For the 
continuation techniques used by software tools like AUTO [39], 
however, the inverse Jacobian plays a central role, so they cannot 
be applied to the system unless this degeneracy is resolved. 
Furthermore the phase space structure of such nonh^pcrbolic 
systems can be changed with arbitrarily small perturbations which 
is why they are called structurally unstable [40]. Note that the 
Unear dependence can be avoided when the rate equation for V 
contains an additional current with a fixed reversal potential 
breaking the symmetry. Such strictly speaking unphysical currents 
are indeed often included in neuronal ion-based models 
[17,24,33,34], but we wiU rather make use of the symmetry and 
eliminate one linearly dependent variable. 

The physiological view on the instability should be as follows. 
Assume that the system is in its physiological equilibrium and then 
apply a constant current lapp to the voltage rate eq. (23). Then eqs. 
(17)-(19) and (23) imply that the equilibrium conditions V = 0 and 
Ki = Noi = Cli = 0 are contradictory, so the equilibrium wiU 
vanish even for arbitrarily small currents. In fact for any constant 
and positi\i' lapp, tlu; system will evolve in a highly non- 
physiological manner with Kj, Nui and C4 slowly tending to zero. 

To avoid even the theoretical possibility of such behaviour we 
will now use eq. (28) to reduce the system and thereby make it 
structurally stable. We can, for example, eliminate V and express 
it in terms of the ICS ion concentrations rather than treating V as 
an independent dynamical variable: 



^(v--^{Nai+Ki-Cli)\=0 

dt \ CmJ 



-Kf^-Ch + Clf') 



This was also done in Ref [l.^i]. The physiological meaning of this 
reduction is simply that the possibility of unspecified applied 
currents is ruled out. For instance, a perturbation on the sodium 



rate eq. (17) should be interpreted as a sodium current. The above 
constraint describes the simultaneous effect on V. It would be 
equivalent to apply perturbations to eq. (17) and eq. (23) 
consistently to model the full effect of an applied sodium current, 
so the additional constraint should be seen as a consistency 
condition. (The curves in Fig. 1(a) were computed for a sodium 
current pulse.) This consistency rule does not at all change the 
dynamics unless unspecified currents are applied, and even then it 
practically does not change the dynamics, because any deviation in 
ion concentrations scales with y and is hence negligible. The 
structural instability is thus a rather formal feature of the 
degenerate model and we remark that its physiological equilibrium 
is nevertheless stationary. Instabilities that lead to an unphysio- 
logical drift of ion concentrations for very long simulation times 
have been reported and resolved in cardiac cell models [7,8]. Our 
case is different though, because the physiological state is a 
stationary one and the response to moderate stimulation is 
physiologically realistic. 

For the bifurcation analyses presented in this paper we have 
eliminated Noi rather than V for numerical reasons. This is 
completely equivalent, because we only vary the pump rate and 
morphological parameters. So in our reduction we have replaced 
rate eq. (1 7) by the following constraint: 



Nai = Na 



,{0) 



+ 



{V- 



■Cl] 



<o) 



(29) 



The model is then defined by the rate eqs. (2), (18), (19) and (23) 
and die constraint eqs. (15), (16), (24)-(26) and (29). 

Bifurcation analysis. We have used the continuation tool 
AUTO [39] to follow the polarized fixed point of the system under 
variation of the maximal pump rate p. Stability changes and the 
creation of stable or unstable limit cycles are detected by the 
software which helps us to interpret the dynamical behaviour. For 
a better overview we will extend our bifurcation analysis even 
beyond the physiologically relevant range. The fuU bifurcation 
diagram is presented in Fig. 2. 

In the (p, K)-plane the fixed point continuation yields a smooth 
z— shaped cur\'e where unstable sections are dashed. The 
physiological equilibrium is marked by a green square. For higher 
pump rates the equilibrium remains stable and becomes slightiy 
hyperpolarized. If p is decreased the physiological equilibrium 
collides with a saddle point at p^^p^ = 0.894006 (xA/cm-^ in a 
saddle-node bifurcation (limit point, LP). In a LP the stability of a 
fixed point chaiigcs in one direction (zero— eigenvalue bifurcation). 
Thus after LPl the fixed point is a saddle point with one unstable 
direction. In a Hopf bifurcation (HB) at pjjgj =29.2336 (tA/cm^ 
two more directions become unstable. Via another LP at 
Pz,P2 = 34.5299 |iA/cm^ the last stable direction switches to 
unstable and the saddle becomes an unstable node. In HBs at 
Pirn = 33.7285 |iA/cm^ and p;y^3 =24.6269 nA/cm^ the fixed 
point becomes a saddle and a stable depolarized focus, respec- 
tively. The stability is indicated by the (n_,«+)-tuples along the 
fixed point curve with denoting the number of stable and 

unstable directions. 

In every HB a limit cycle is created. Our model only contains 
unstable limit cycles that are created in subcritical HBs. In the 
diagram they are represented by their extremal V values. Such 
unstable limit cycles are not directiy observable, but in the bistable 
regime they can play a role for the threshold behaviour for the 
transition from one fixed point to the other. All limit cycles in the 
model disappear in homochnic bifurcations (HOM). In a HOM a 
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Figure 2. Bifurcation diagram of tlie ion-based model. Bifurcations are marked by red circles, the physiological equilibrium by a green square. 
Following the z-shaped fixed point characteristic from below there are two saddle-node bifurcations (limit point, LP) at p = 0.894006 |iA/cm^ and 
p = 34.5299 tiA/cm^ and three subcritical Hopf bifurcations (HB) at p = 29.2336 |xA/cm^, p = 33.7285 nA/cm^ and p = 24.6269 nA/cm^. The limit 
cycles created in HB1, HB2 and HB3 disappear in homoclinic bifurcations (HOiVl) at p = 27.6463 |iA/cm^, p = 33.7027 |iA/cm^ and 
p = 1.024291 |xA/cm^, respectively. The second LP and the second HB together with the HOIVI of limit cycles occur in a very narrow parameter 
range (see blow-up inset). The number of stable («_) and unstable (« + ) directions of the fixed point is indicated by the )-tuples. There is 

bistability of a physiological state and a depolarized state with largely reduced ion concentration gradients between p = 8.94006 |iA/cm^ and 
p = 24.6269 nA/cm^. 
doi:1 0.1 371 /journal.pcbi.1 003551 .g002 



limit cycle collides with a .saddle. When it touches the saddle it 
becomes a liomoclinic cycle of iiifmite period. After the bifurcation 
the limit cycle does not exist any more. The limit cycles created in 
HBl, HB2 and HB3 disappear in HOMs at p = 27.6463 nA/cm^, 
p = 33.7027 nA/cm^ and p= 1.024291 [lA/cm^. Tlie limit cycle 
emanating from HBl collides with the upper (i.e., less polarized) 
saddle, for the other two HOMs the situation is clear, because 
there is only one saddle available. Since the limit cycles are all 
unstable these bifurcation details are physiologically irrelevant, but 
mentioned for completeness. 

This bifurcation analysis shows that our model is bistable for a 
large range of pump rates p^pj <p< Phb3 ■ Strongly depolarized 
and electrically inactive states of neurons with nearly vanishing ion 
concentration gradients have been reported in pathological states 
[11,13], but in real systems such free energy-starvation (FES) is 
not stable. In the below section we show how this bistability can be 
resolved. 

Ionic excitability. We wiU now briefly show how the above 
analyzed model can be modified such that the unphysiological 
bistability turns into excitability of ion dynamics. For this we follow 
[24,33] and include an additional regulation term for extracellular 
potassium. This means that becomes an independent 

dynamical variable and the constraint eq. (25) must be replaced 
by its rate equation: 

^ = + ^^^^ 

The regulation term Ireg can be interpreted as a diffusive coupUng 
to an extracellular potassium bath or as a phenomenological 
buffering term. It takes the following form: 



Ireg = KKreg — Ke), (31) 

where K^g is the potassium concentration of an infinite bath 
reservoir coupled to the neuron or a characteristic parameter for 
glial buffering, and /I is a rate constant (values given in Table 3). 
Kri,g takes values of physiological potassium concentrations and 
hence stabilizes the physiological equilibrium. This is how lyeg 
regulates ion homeostasis and destabilizes the energy-starved 
state. 

If we now stimulate the system with a current pulse or 
temporarily switch off the pump as we did in Fig. 1 the system no 
longer remains in the depolarized state, but repolarizes after a long 
transient state of FES (see Fig. 3). After the repolarization ion 
concentrations start to recover from FES. FuU recovery to the 
initial physiological values is an asymptotic process which takes 
very long (about two hours), but the neuron is back to normal 
functioning already after nine to ten minutes. Similar dynamics is 
described in numerical [17,34] and experimental SD models. 

The bistable and excitable dynamics can be nicely compared in 
a projection of the respective trajectories onto the (iVa, ,.fir,)-plane. 



Table 3. Buffering parameters. 







Name 


Value & unit 


Description 


A 


2.7e-5/msec 


regulation rate 


l^reg 


4 mMol// 


regulation level 



doi:l 0.1 371 /journal.pcbi.1 003551 .t003 
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Figure 3. Upper panels: Membrane and Nernst potentials, lower panel: Ion concentrations vs time, (a) Stimulation with the same, but earlier 
applied, current pulse as in Fig. 1 (a). Due to the additional potassium regulation the system returns to the physiological equilibrium after an approximately 
60 sec lasting FES and subsequent hyperpolarization. (b) Similar dynamics as in (a) is observed for a temporary pump switch-off like in Fig. 1(b). 
doi:1 0.1 371/journal.pcbi.1 003551 .g003 



For the bistable model the conditions K, = 0 and iVfl, = 0 define 
three-dimensional hypersurfaces called nullclines. Adding the 
necessary fixed point conditions on the remaining dynamical 
variables, namely V = Eci and n = nao{V), allows us to specify 
curves that represent these nullclines and only depend on Na, and 
Ki. In the buSered model is another dynamical variable and its 
frxed point condition is Ke = Kreg. Electroneutrality and mass 
conservation imply that certain (A'fl,-,A^,)-combinations would lead 
to a negative C/, or K^. In the plot, these unphysiological 
configurations are shaded. 

In Fig. 4, we see that the bistable model has three nuUcline 
intersections, i.e., fixed points, while the buffering term deforms 
the nullclines so that only one stable fixed point remains. In the 
bistable case an initial current pulse stimulation (dashed part of 
trajectory) drives the system into the basin of attraction of the 
FES-state, which it then asymptotically approaches (solid line). 
After the same stimulation the buffered system performs a large 
excursion in phase space with extremal ion concentrations 
comparable to FES, but eventually returns to the physiological 
equilibrium. This large excursion in the ionic variables character- 
izes what we refer to as ionic excitability or excitability of ion 
homeostasis. The simulations presented in this section support the 
hypothesis that it is caused by the bistabUity of the unbuffered 
model. Note that intersections of nuUcline curves and trajectories 
do not have to be horizontal or vertical since they may (and do) 
differ in the non-ionic variables. The main purpose of the nuUcline 
curves is to indicate the existence and location of fixed points. 



Robustness of results 

The ion-based model we have analysed so far has been 
motivated as a natural extension of the Hodgkin-Huxley 
membrane model. However, there are different variants of ion- 
based models [17-24,32] that use different pump and current 
models, ion content, and ion channels. We wUl hence address the 
question how general our results are in this respect. Furthermore 
we vary the geometry-dependent parameters (membrane surface 
and extraceUular volume fraction) continuously to test their effect 
on the phase space, too. 

Model variants. As we noted before, transmembrane cur- 
rents are more accurately described as Goldman-Hodgkin-Katz 
(GHK) rather than Nernst currents, even though we prefer the 
latter. It is hence important to check which difference the choice of 
current model makes. To generalize the Nernst currents in eqs. 
(11)-(13) to GHK currents we assume that both models have the 
same steady state currents under physiological equUibrium 
conditions. The GHK version of the sodium current is 

Noi-Na, exp( - K/26.64) ^^^^ 
l-exp(-K/26.64) 

with membrane permeabilities and P^,, instead of conduc- 
tances. To compute these permeabUities we set the GHK current 
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Figure 4. Projection of the trajectories corresponding to Fig. 1a and Fig. 3a. The shaded regions indicate unphysiological {Nai,Kf)- 
combinations that imply negative CI, (lower left region) or negative K^. (upper region in left plot). Stable and unstable fixed points are marked by 
solid and open circles, (a) In the bistable case an initial stimulation (dashed line) leads to large subsequent changes in ion concentrations that 
terminate in the second fixed point of the system, (b) The excitable motion starts very similar to case (a), but after reaching the extremal 
concentration values the system slowly returns to its initial state. 
doi:1 0.1 371/journal.pcbi.1 003551 .g004 



equal to its NeriLstian counterpart 

lN.+ =(g'Na+g'Na'^'h)iV-E^,) (33) 

for the equilibrium conditions given in Table 1. This leads to a 
common conversion factor from conductances g"^^ to permeabil- 
ities Pji'^. With this ansatz we obtain conversion factors for the 
three different ion species that lead to the conductances listed in 
Table 4. 

There is also a certain freedom in the choice of a pump model. 
It is a general feature of Na"'"/K"'" pumps that their activity is 
enhanced by the elevation of ECS potassium and ICS sodium. Still 
different models exist, and to investigate the role of the particular 
pump model we replace the pump from eq. (21), now referred to 
as Ip^A, with the following one from [17]: 

Ips(Na>,K,) = Ps(^l + ^^ ^^^^ 

In order to retain the equilibrium at K = — 68 mV we have to set 
the maximum pump current to = 5.72 |iA/cm^. This is slightly 
higher than the previous pump value (p^ = 5.25 |iA/cm^), but in 
the same range. 

From the rate equations of the HH membrane model (see Sec. 
Model) it is obvious that the chloride leak current stabilizes the 
equilibrium membrane potential. To test its stabilizing effect in the 



Table 4. Mennbrane pernneabilities for GHK current. 





Name 


Value & unit 


Description 


P' 


0.0264 |im/sec 


leak sodium permeability 




150.77 |im/sec 


gated sodium permeability 




0.0169 |im/sec 


leak potassium permeability 


P^ 


13.488 nm/sec 


gated potassium permeability 


P' 


0.0521 nm/sec 


leak chloride permeability 



doi:l 0.1 371 /journal.pcbi.l 003551 .t004 



context of ion-based modeling we compare models that either do 
or do not contain this current. We are further interested in the 
question whether membrane excitability and ion bistability are 
related. Therefore also the effect of in- and excluding active ion 
channels is tested. 

In this section we wUl only discuss fixed points and their 
stability, but not the unstable limit cycles belonging to HBs. In 
Fig. 5 the fixed point continuation curves for aU combinations of 
current model (Neriist or GHK), pump choice {Ip^A or Ip s) and the 
respective in- and exclusion of chloride and active ion channels 
channel are shown. Each panel (a)-(d) contains all continuation 
curves for a given choice of pump and current model. For those 
models that are bistable for certain pump rates an overview of the 
different dynamical regimes is presented in Fig. 6. It shows the 
parameter ranges for bistability and for monostability of a 
physiological state or FES. 

The most striking result of this bifurcation analysis is that this 
bistability occurs in all models with gated ion channels, but not in 
any model with only leak channels (grey-shaded graphs in the 
insets of Fig. 5). The comparison of any model with active gates 
and its leak-only counterpart shows that whenever the physiolog- 
ical equilibrium of the first one exists it is identical to the 
equilibrium of the latter one. While the physiological state 
disappears in a LP for all bistable models at small pumping, the 
fixed points of the leak models remain stable, but depolarize 
drastically for further decreasing pump rates until the Donnan 
equilibrium for p^ ^ = 0 is reached. The absence of the second 
fixed point in leak-only models is plausible if we consider Fig. 1 
again. The depolarized state is characterized by large ion 
concentrations Ke and Noj which implies an increased pump 
current (see eq. (21)). Since the differences between the Nernst 
potentials and the membrane potential are even smaller in the 
depolarized state, higher, hence gated, conductances are required 
to compensate for the pump currents and maintain the 
depolarized state. Besides the requirement of active ion channels, 
the bistability is a very robust feature of these simple ion-based 
models. 

Let us now consider the effect of the different model features on 
the minimal physiological pump rate, i.e., the pump rate required 
for a stable physiological fixed point, and the recovery pump rate 
that destabilizes the depolarized state of FES and allows the 
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Figure 5. Bifurcation diagrams of fixed points for different models. The effects of chloride and active ion channels are compared for each of 
the four possible pump (A vs B) and current model (Nernst vs GHK) combinations. The physiological equilibrium for normal pump rates 
{pj = 5.25 |xA/cni^ and pg = 5.72 |iA/cm^) is marked by a green square. The model from fig. 2 is marked by a star. The value is the same with and 
without chloride or active channels. Insets show the bifurcation diagrams for low pump rates (p^ g <5 |xA/cni-). Fixed point lines for models without 
active ion channels are shaded (see insets). Note the different scales on the main figures, insets are for the same range in each panel. 
doi:10.1371/journal.pcbi.1003551.g005 
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Figure 6. Overview of tlie parameter regimes for bistabiiity, polarized and depolarized stability for different models (1-8). The 

change from the monostable depolarized regime to bistability (red to orange) defines the minimal physiological pump rate, i.e., the pump rate 
required for the existence of a polarized fixed point. The line separating the bistable from the monostable polarized regime (orange to green) defines 
the minimal recovery pump rate, i.e., the pump rate required to return from the depolarized fixed point to the polarized equilibrium. The model from 
Sec. Model is marked by a star. 
doi:1 0.1 371 /journal.pcbi.1 003551 .g006 



neuron to return to physiological conditions. These rates are the 
lower and upper limit of the bistable regime, and low values are 
physiologically desirable. 

In Fig. 6 we see that pump model A, GHK currents and 
chloride each lead to a lower minimal physiological (see the inset) 
and a lower recovery pump rate than pump model B, Nernst 
currents and the exclusion of chloride. Quantitative diflFerences 
should be noted, though. The inset of Fig. 6 shows that all models 
with pump A have lower minimal physiological pump rates than 
models with pump B. So the stability of the physiological 
equilibrium with respect to pump strength reduction depends 
mostly on the choice of pump. On the other hand four of the five 
lowest recovery pump rates are from models that include chloride. 
In fact, it is only the combination of both, the GHK current model 
and pump A, that make the recovery threshold of the chloride 
excluding model 7 slightly lower than that of the chloride 
including model 2. However, one should note that even the 
lowest recovery pump rate is as high as p= 14.3 |iA/cm^ (model 
8). This is stiU an almost threefold increase of the normal rate. So 
even if we assume pump enhancement due to additional 
mechanisms, for example increased cerebral blood flow, the 
threshold for recovery from FES seems to be too high. Thus it is 
true for a large class of ion — based neuron models that realistic 
neuronal homeostasis cannot rely on Na''"/K"*"-ATPase alone, but 
rather on a combination of ion pumps and further regulation 
mechanisms like glial buffering. 

There is another effect of chloride to be pointed out. In Fig. 5 
we see that it raises the Donnan equilibrium potential (see 
potentials at p = OpA/cm^) significantiy. To understand this effect 
note that without chloride electroneutrality forces the sum of AK^ 
and ANOe to be zero, while the presence of the decreasing anion 
species C4 implies A{K^, + Nae)<0. According to eq. (14) this 
leads to lower Donnan equilibrium Nernst potentials Ek and E^a, 



and consequently to a lower membrane potential. Since the 
conditions of FES for physiological pump rate values are very close 
to the Donnan equilibrium, this depolarized fixed point is shifted 
in the same way. 

The effect of the current model on the two characteristic pump 
rates is less pronounced than that of chloride or the pump choice. 
It lowers the minimal physiological rate more than chloride, but 
not as much as a pump change from model B to A. Its effect on the 
recovery threshold is the weakest. 

While above we describe and investigate minimal Hodgkin- 
Huxley model variants to obtain SD behavior in the simplest 
neuron model types, in the current literature biophysically much 
more detailed neuron models have been developed for this 
phenomenon. We do not intend to investigate such detailed 
models thoroughly, but as an example that further demonstrates 
the robustness of our results, we also replicate our results with a 
much more detailed membrane model as first described by Kager 
et al. [17]. This detailed model contains five different gated ion 
channels (transient and persistent sodium, delayed rectifier and 
transient potassium, and NMDA receptor gated currents) and has 
been used intensively to study spreading depolarizations and 
seizure-like activity. In fact, one modification is required so that 
we can replicate our previous results. The detailed model contains 
an unphysiological so-called 'fixed leak' current 

Ileak/=gleak/iV + 70) (35) 

that has a fixed reversal potential of —70 mV and no associated 
ion species. This current only enters the rate equation for the 
membrane potential V and thereby implies that K = — 70 mV is a 
necessary fixed point condition. In other words, the type of 
depolarized fixed point that we have found in the simpler model is 
ruled out by this current. If we, however, replace this unphysical 
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Figure 7. Bifurcation diagram of the fixed points for the Kager et ai. model [17]. The physiological equilibrium is at p= 13 |AA/cm^, the 

minimal physiological pump rate is p = 9.8 |jA/cm-, and the recovery rate is p = 107 |iA/cm^. The limit cycle emanating from the HB undergoes four 
saddle-node bifurcations of limit cycles (indicated by the stability changes, but not explicitly labled) before it disappears in a homoclinic bifurcation 
(HOIVl). Like in Fig. 5 the fixed point line for the corresponding leak-only model is shaded. Its value at p = 0 |j.A/cm^ indicates the Donnan 
equilibrium. 

doi:10.1371/journal.pcbi.1003551.g007 



current with a chloride leak current as in our model (see eqs. (13), 
(19), (23)) and furthermore neglect the glial bufiering model use in 
Ref [17], we find the same type of bistability as in our model. 

The fixed point continuation of this model (for a complete list of 
rate equations see [17,34]) in Fig. 7 shows that again FES 
conditions and a physiological state coexist for a large range of 
pump rates. This model has slightiy different leak conductances 
and equilibrium ion concentrations, and consequently the 
characteristic pump rates also differ from ours. However, the 
only important thing to note here is that the recovery pump rate 
defined by the subcritical Hopf bifurcation of the upper fixed point 
branch is large compared to the physiological value (marked by 
the green square), so also in this rather different membrane model 
recovery from FES due to pump enhancement is practically 
impossible. The limit cycle continuation also bears a strong 
similarity to the one in Fig. 2 with the only main difierence being 
that the limit cycles of the detailed model are stable in two narrow 
parameter regimes (see solid circles in Fig. 7). We remark that the 
physiologically irrelevant unstable frxed point branches in Fig. 7 do 
not connect in a saddle-node bifurcation, but saturate for very 
high pump rates. The occurrence of the same type of bistability in 
a Hodgkin-Huxley-based ion model and this very detailed one, 
and also the similarity of the SD trajectories in Fig. 3 to those 
presented in [17], support the physiological relevance of the 
minimal ion-based ansatz that we developed and follow in this 
paper. Moreover, we assume that bistability is an universal feature 
also for other more detailed membrane model, which are yet more 
elaborate variants of the model described by Kager et al. [17]. We 
wiU hence use the model from Sec. Model for further investiga- 
tions. 

Variation of membrane surface and extracellular volume 
fraction. After the overview of different variants of ion content, 
ion channels, pumps and current models we finally address the 
role of the neuron geometry. Therefore we vary the membrane 



surface and the extracellular volume fraction in the model from 
Sec. Model. For the surface variation we introduce the relative 
surface size parameter Xa '^'^d replace A,„ with A„,x^ which 
implies the replacement (see eq. (20)) 

y^yxA 

wherever y occurs, i.e., in the ion rate eqs. (18) and (19) and the 
sodium constraint eq. (29). The extracellular volume fraction, 
typically denoted as /, is defined as 

/:=— => oje=/co,o, and co, = (1 

where Wjoi = co/ + (t>e is the total volume of the system. When / is 
varied, the above expressions for £0,/^ must be inserted in both ion 
rate eqs. (18) and (19), and in all ion constraint eq. (24)-(26) and 
(29). The surface parameter is varied from 0.1 to 10,/ is varied 
from 2% to 50%. The standard values of these parameters are 
Xa = ^ and / = 25% and parameters are understood to take these 
values when they are not varied. We start from the bifurcation 
diagram of Fig. 2 and perform two-parameter continuations of the 
detected bifurcations to find out how the membrane surface and 
the volume ratio change the bistable regime. The (p,/^)- and 
(pj^) — continuation curves are shown in the left and right plot of 
Fig. 8. 

We see that the Xa variation has hardly any effect on the 
bifurcation values of p. This can be understood from the structure 
of the model. The fixed point curve is defined by setting the rate 
eqs. (2), (18), (19), (23) to zero and the constraint eqs. (15), (16), 
(24)-(26) and (29). When Xa is varied the only modification to 
these conditions is in eq. (29) (sodium constraint). But this 
modification is of order 0(10^') and does practically not affect 
the shape of the frxed point curve, so the limit point bifurcations 
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Figure 8. Two-parameter continuations of the fixed point bifurcations of Fig. 2. In the left plot the dimensionless surface size parameter 
is varied, in the right plot the extracellular volume fraction / is changed. The insets show the LP1 curves that mark the minimal physiological pump 
rate. The pump rates for which the system is bistable range from the LP1 to the HB3. The HB3 pump rate is required to repolarize a neuron that is in 
the depolarized equilibrium. The parameter (left plot plus inset) almost does not change the stability of the system, but / (right plot) reduces the 
recovery pump rate significantly. The inset shows that the minimal physiological pump rate is much less affected. In each plot and inset the standard 
parameter value is indicated by the light-blue vertical line. 
doi:1 0.1 371 /journal.pcbi.1 003551 .g008 



LP 1 and LP2 are almost not changed. Hopf bifurcation could be 
shifted, but a rescaling of the (initial and dynamical) ion 
concentrations by transforms the rate and constraint equations 
such that Xa o"ly appears in the pump currents. Their derivatives 
are then multiplied by Xa^ but for all HBs the pumps are saturated 
and hence Xa does not contribute to the Jacobian. The variation of 
/, however, does change the width (with respect to p) and the 
threshold values of the bistable regime. A small value of / 
(corresponding to a small extracellular space) reduces the recovery 
pump rate, and also increases the minimal physiological pump 
rate. This means that both, depolarization and recovery, are 
enhanced. However, the minimal physiological pump rate is much 
less affected than the recovery pump rate, so basically a big cell 
volume supports recovery from the depolarized state. It is known 
that in spreading depression (SD), where metastable depolarized 
states that resemble the energy-starved fixed point of our model 
occur, the osmotic imbalance of ICS and ECS ion concentrations 
leads to a water influx that makes the cells swell. Our analysis 
shows that such a process helps the neuron to return to its 
physiological equilibrium. Extracellular volume fractions of down 
to 4% are reported in SD, but even for such extreme volume 
fractions the required recovery pump rate is too high for pump 
driven recovery of the neuron (see the lowest value for PnBi ™ the 
right plot of Fig. 8). We remark that the bifurcation curves in Fig. 8 
do not saturate for />50%, but, except from the LPl-curve that 
remains very low, bend down, probably due to an approximate 
symmetry of ICS and ECS ion concentration dynamics. In 
summary also the analysis of different cell geometries confirms that 
ion homeostasis cannot be provided by Na"'"/K"'" pumps alone. For 
example, in computational models of dynamically changing pump 
rates due to oxygen consumption maximal rates of twice the 
physiological values are considered [22]. 

Discussion 

Computational neuroscience complements experimental and 
clinical neuroscience. Simulations help to interpret data and guide 
a principal understanding of the nervous systems in both health 
and disease. The HH-formulation of excitability was "so 
spectacularly successful that, paradoxically, it created an unreal- 
istic expectation for its rapid application elsewhere" as Noble 



remarked [9] . While his statement refers to modeling of cardiac 
cells it certainly holds true also for neurological diseases and brain 
injury [11,13]. In both fields, the incorporation of the Na"'"/K"'' 
pump in the original excitability paradigm formulated by Hodgkin 
and Huxley is of major importance. The fundamental structure of 
such models has to our knowledge not been exploited in 
neuroscience beyond merely modulating spiking in epileptiform 
activity [23,24] or in models that have energy-starved states [17- 
22,32] yet without investigating the fundamental bifurcation 
structure. 

As we stressed in the introduction, this extension of the original 
HH model enforces a physical or rather thermodynamical 
perspective, which was, of course, the starting point of Hodgkin 
and Huxley, too. For instance, we also considered the Goldman— 
Hodgkin-Katz (GHK) current equation which is derived from the 
constant field assumption applied to the Nernst-Planck equation 
of electrodifFusion. Electroneutrality is important to consider, as 
can be seen by the indirect insertion of impermeable counter 
anions only reflected by observing a thermodynamic Donnan 
equilibrium. Furthermore, a thermodynamic description of 
osmotic pressure (which would require a direct insertion of a 
concentration A"^ of a counter anion with valence n) and 
corresponding changes in cell volume can be included. 

There are further physical mechanisms that may alter the 
dynamics in biophysical ion-based models. At the same time, we 
have to avoid "an excruciating abundance of detail in some 
aspects, whilst other important facets [...] can only be guessed" 
[41], like using various new currents but guessing the correct value 
of the valence n of an impermeable counter anion. For this reason, 
we decided to use the original ion currents from the HH model. 
The comparison of our results to a physilogically more realistic 
and much more detailed membrane model in Fig. 7 support the 
assumption that the basic structure will not be changed by just 
adding or modifying gating. This question has also been addressed 
experimentally and in simulations by showing that only the 
simultaneous blockade of all known major cation inward currents 
did prevent hypoxia-induced depolarization with FES [17,42]. In 
the model [17], five different Na^ currents were investigated. Of 
course, to apply our model to a particular pathological condition, 
like migraine which is a channelopathy [43,44] (disease caused by 
modified gating), these details will become important. This can 
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easily be incorporated in future investigations. Moreover, note that 
changes in cell volume, which are very important in brain injuries, 
are in this study only treated by varying it as a parameter. 

Our bifurcation analysis shows that a whole class of minimal 
ion-based models is bistable for a large range of pump rates 
(Plpi < P <Phb3)- Bistable dynamics was suggested by Hodgkin to 
explain spreading depression [11,12], and a corresponding model 
has been investigated mathematically by Huxley but never been 
published (c£ Ref [45]). Dahlem and Miiller suggested to extend 
this ad hoc approach, i.e., a single so-called activator variable with 
a bistable cubic rate function, by including an inhibitory 
mechanism in form of an inhibitor species with a linear rate 
function coupled to the activator [45] . This, of course, leads to the 
well known FitzHugh-Nagumo paradigm of excitability type II 
[46,47], that is, excitability caused by a Hopf bifurcation [48], but 
should not be mistaken as a modification of conductance-based 
excitability in form of HH-type model in the 'first generation' and 
the interpretation as an equivalent electrical circuit. FitzHugh used 
his equation in this way, he investigated a long plateau as seen in 
cardiac action potentials. Dahlem and MtiEer suggested to use the 
same mathematical structure of an activator-inhibitor type model 
[45] to describe a fundamental new physiological mechanism of 
ionic excitability that originates from bistable ion dynamics. Our 
current results provide the missing link between this ad hoc 
activator-inhibitor approach, which has been widely used in 
migraine and stroke pathophysiology [45,49-58], and biophysical 
plausible models. The major result from this link is the new 
interpretation of the physiological origin of the proposed inhibitor\' 
variable [45]. We wrongly interpreted it as being related to the 
pump rate [49,51,53,57]. 

As our ion-based model shows bistable dynamics, we see it as 
essentially capturing the activator dynamics of an excitable system 
and briefly show in Figs. 2 and 3 that it can be transformed into 
such a system by the introduction of a inhibitory process. Vice 
versa, excitable systems can be reduced to bistable dynamics by 
singular perturbation methods. Such reductions are referred to as 
a threshold reduction. From this perspective our model can be 
interpreted as the threshold reduction of an excitable system, and 
we conclude that without contact to an ion bath, physically 
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